! yasm -a x86 -f elf dec.s -m amd64 &&  gfortran45 dec.o test_dec.f90 -o ~/bin/test_dec
program list2hist
  implicit none
  interface 
     subroutine decs(e,xyz,etc)
       integer*8,intent(in)::e
       integer*4,intent(out)::xyz(3)
       integer*4,intent(out)::etc(4)
     end subroutine decs
  end interface      
  integer,parameter::nx=256
  integer,parameter::ny=256
  integer,parameter::nxb=2
  integer,parameter::nyb=2
!!$  integer,parameter::nx=512
!!$  integer,parameter::ny=512
!!$  integer,parameter::nxb=1
!!$  integer,parameter::nyb=1
  integer,parameter::nf=5
  !
  integer*4,target::xyz(3),etc(4)

  integer*4,allocatable::buf4(:)
  character*256 fin,title
  integer*4 istat,len
  integer*8 n,ie,b1,b2,nb
  integer*8 fsize
  integer,pointer::ix,iy
  integer jx,jy
  !
  logical lo_found
  integer nok
  integer*8 ev
  pointer(pev,ev)
  integer nlost,nnull
  integer hist(nx,ny)
  integer if
  integer istep
  integer pmtid
  !
  title=""
  n=command_argument_count()
  if(n<2) stop "gimme pmtid flist [title]"
  call get_command_argument(1,fin,len)
  read(fin,*,iostat=istat) pmtid
  if(istat/=0) stop "*** invalid pmtid"
  call get_command_argument(2,fin,len)
  fin=trim(adjustl(fin))
  if(n>=3) call get_command_argument(3,title,len)
  title=adjustl(title)
  !
  open(unit=1,file=fin,status="old",iostat=istat)
  if(istat/=0) stop "*** cannot open flist"

  ix=>xyz(1)
  iy=>xyz(2)

  if=0
  hist=0
  istep=0
  file: do
     read(1,"(a)",iostat=istat) fin
     if(istat/=0) exit
     open(unit=2,file=fin,status="old",access="stream",iostat=istat)
     if(istat/=0) then
        write(*,*) "cannot open file: "//trim(fin)
        stop
     end if
     if=if+1
     call fseek(2,0,2) ! seek end
     call ftell(2,fsize)
     call fseek(2,0,0) ! seek set
     n=fsize/4
     if(allocated(buf4)) deallocate(buf4)
     allocate(buf4(n),stat=istat)
     if(istat/=0) stop "*** cannot allocate buf"
     b1=1
     nb=1024*1024*1024/16
     
     list: do ie=1,n/nb+1
        b2=min(b1+nb,n)
        read(2,iostat=istat) buf4(b1:b2)    
        b1=b2+1
        if(istat/=0) then
           write(*,*) ie,istat
           stop "*** error during read."
        end if
     end do list
     close(2)
     
     lo_found=.false.
     nok=0
     nlost=0
     nnull=0
     event: do ie=1,n
        if(buf4(ie)==0) then
           nnull=nnull+1
           cycle
        end if
        if(lo_found) then
           if(get_flg(buf4(ie))==0) then
              ! pair found
              nok=nok+1
              lo_found=.false.
              pev=loc(buf4(ie-1))
              call decs(ev,xyz,etc)
              if(pmtid>0.and.pmtid/=etc(2)+1) cycle
!              write(*,200) xyz(1),xyz(2),xyz(3),etc(1),etc(2),etc(3),etc(4)
              jy=dble(iy)/dble(nyb)+1
              jx=dble(ix)/dble(nxb)+1
              if(jx>0.and.jx<=nx.and.jy>0.and.jy<=ny)&
                   hist(jx,jy)=hist(jx,jy)+1
           ELSE
              nlost=nlost+1
              lo_found=.true.
           end if
        else        
           if(get_flg(buf4(ie))==1) then
              lo_found=.true.
           else
              nlost=nlost+1
           end if
        end if
     end do event
     if(if==nf) then
        istep=istep+1
        call write_hist
        hist=0
        if=0
     end if
  end do file
  
contains

  integer function get_flg(e)
    integer*4,intent(in)::e
    get_flg=ishft(e,-30)
  end function get_flg
  
  subroutine write_header
    character*512 hdr
    write(hdr,*) istep*nf
    hdr="Processed="//trim(adjustl(hdr))
    write(3,100) trim(title),trim(hdr),nx,nx,ny,ny,maxval(hist)
    100 format("$ DATA=CONTOUR",/,&
             "%toplabel=",a,/,&
             "%subtitle=",a,/,&
             "%contstyle=2 nsteps=50",/,&
             "%XMIN=1",x,"XMAX=",i4,x,"NX=",i4,/,&
             "%YMIN=1",x,"YMAX=",i4,x,"NY=",i4,/,&
             "%cmin=0 cmax=",i9)
  end subroutine write_header

  subroutine write_hist
    character*512 fout
    write(fout,*) istep
    fout=trim(adjustl(fout))//".mtv"
    open(unit=3,file=fout,iostat=istat,status="replace")
    if(istat/=0) stop "*** cannot open fout"
    call write_header
    do jy=1,ny
       do jx=1,nx
          write(3,*) hist(jx,jy)
       end do
    end do
    close(3)
  end subroutine write_hist

end program list2hist
